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In  reactor  systems,  the  reaction  rates  are  calculated  from  a velocity  space  integral  over  the 
reaction  cross-section  times  the  distribution-functions  of  the  reacting  species.  In  a plasma,  the 
shape  of  the  distribution  functions  themselves  is  determined  by  solving  a diffusion  plus  convection 
problem  in  velocity  space.  This  review  describes  briefly  the  physical  processes  of  central  interest 
to  such  a description  and  the  mathematical  formulation  of  the  problem.  It  presents  the 
numerical  methods  which  have  been  u.sed  in  such  calculations  by  various  authors.  Optimization 
on  a vector  computer  (the  Texas  Instruments  ASC)  is  described.  Finally,  .some  indication  is 
given  of  what  may  be  expected  as  reactor  systems  are  treated  in  more  detail.  ^ 
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Velocity-Space  Methods  for  Reactor  Plasmas 


INTRODl  t TION 

The  Fokker-Planek  Tguation  is  one  of  a number  of  specialized  tools  that  can 
give  you  a handle  on  the  dynamics  of  reactor  plasmas.  Like  any  tool,  it  is  more  useful  m 
some  cases  than  m otners  fhe  emphasis  here  is  on  providing  the  potential  user  with  the 
onent.ition  he  needs  to  determine  the  usefulness  of  these  methods  for  him. 

This  paper  presents  a brief  review  of  the  analytic  and  numerical  development  of 
the  I okker-Planck  equation,  and  an  introduction  to  recent  work,  much  of  which  is 
presently  unpublished.  SI  units  are  used  throughout  The  presentation  is  reasonably  self- 
contained.  but  no  attempt  is  made  to  present  the  details  of  derivations  or  applications. 
References  are  provided  lor  those  who  wish  to  examine  these  details.  The  succeeding  sec- 
tions are  organized  by  discipline,  so  that  the  reader  with  particular  interests  may  quickly 
linil  that  part  ol  the  material  which  interests  him  most  The  section  on  physics  also  pro- 
vides motivation  for  the  present  development  of  these  methods,  and  the  reactor  section  indi- 
cates the  neeits  for  future  development. 

I he  physics  section  emphasizes  the  difference  between  the  gaseous  and  plasma  state 
It  explains  why  velocity-space  methods  are  necessary  for  plasma  reactors.  It  depicts  several 
practical  physical  problems,  and  explains  how  they  are  represented  in  the  Tokker-Planck 
i quation 

The  mathematics  section  shows  how  the  Fokker-Planck  Lquation  is  derived,  and 
what  assumptions  are  made  It  presents  the  general  form  of  the  equation,  and  described  the 
initial  and  boundarv  values  needed.  Some  analytic  methods  for  solving  it  are  presented 

fhe  numerics  section  describes  several  ways  to  reduce  the  computational  complexity  ol 
the  general  Tokker-Planck  I quation.  Numerical  stability  is  discussed.  Direct  and  transform 
methods  of  solution  are  presented.  The  speed  gains  available  on  vector  computers  are  il- 
lustrated 

I he  applications  section  presents  methods  for  adapting  the  T'okker-Planck  Tguation  to 
experiments  Mirror  and  pinch  systems  are  described.  Tokamak  systems  are  discussed, 
and  the  difficulties  peculiar  to  toroidal  systems  are  explored 

I he  reactor  section  describes  the  diflerences  between  the  present-day  environment  and 
future  reactor  systems  It  presents  some  of  the  problems  associated  with  heating  by  energetic 
reaction  products  It  suggests  some  areas  for  future  development 

P»I>SI(S 

In  reactor  systems,  ihc  nuclear  reaction  rate  depends  on  an  integral 

fi  -JlV  -V'l.r(|V  - V'l  ) /„  (V  ) /^  (V  ) 

whose  main  contributions  are  from  opposite  ends  of  the  distribution  functions 
ol  the  reacting  species  a.h.  This  is  because  the  thermal  velocities  of  a and  b,  and  thus  the 
main  bulk  of  their  distributions,  arc  normally  at  energies  far  below  the  peak  of  the  nuclear 
r reaction  cross-section  Ihc  distribution  functions  t,,  and  arc  rapidly  changing  at  these 

energies,  and  the  integral  K is  very  sensitive  to  such  variations.  Thus  it  is  necessary  to  calcu- 
late and  tf,  quite  accurately 
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This  siiualmn  is  guile  ditTcrenl  from  that  ol  eonvenilonal  lluid  mechanies 
Iherc  only  the  first  levs  momeius  of  the  distribution  functions  (mass,  mtimentum,  energy) 
are  needed,  and  they  may  be  ealeulaled  directly.  The  critical  hypothesis  from  which  fluid 
egualions  are  derived  is  an  assumption  about  the  shape  ol  the  distribution  function,  usual- 
ly that  it  is  nearly  Maxwellian.  This,  in  turn,  depends  on  the  presence  of  some  colli- 
sional  process  that  drives  the  distribution  function  toward  this  slate.  In  conventional  lluids. 
this  assumption  is  quite  reasonable. 


In  the  physics  of  hot  plasmas,  this  assumption  is  not  reasonable.  Collisions  become 
relatively  infrequent  In  fact,  a material  is  considered  to  change  from  its  gaseous  state  to  its 
plasma  state  when  its  plasma  frequency  exceeds  its  collision  frequency,  namely . 
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where  ni,  q.  n and  /'  are  the  mass,  charge,  number  density  and  temperature,  respectively, 
of  the  material,  A is  Holt/mann's  constant,  and  In  \ is  a dimensionless  variable  of  order 
10  [see  Kq.  (11)1  One  cannot  assume  that  the  shape  of  the  distribution  function  is 
known  in  a plasma.  Fortunately,  one  can  assume  that  its  evolution  is  slow,  since  the  col- 
lisions are  infrequent  When  a charged  particle  finds  itself  in  a plasma,  the  other  charged  par- 
ticles around  it  ''shield"  it  Irom  electrostatic  encounters  \n  immediate  consequence  of  this  is 
that  collisional  changes  in  a particle  trajectory  take  place  in  many  small  steps,  rather  than  a 
few  large  ones. 


These  facts  suggest  that  the  evolution  of  the  distribution  is  a dilTusion  process, 
general  form  of  this  diffusion  equation  is 
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where  K and  B are  the  electric  and  magnetic  fields  present  in  the  plasma,  and  K and  L are 
the  diffusion  and  drag  terms.  .Several  specific  properties  of  this  diffusion  process  distinguish  it 
from  more  conventional  cases,  like  heat  flow. 


For  example,  K and  I,  depend  on  all  the  species  present  in  the  plasma,  and  may  be 
anisotropic.  Two  instances  of  this  anisotropy  are  illustrated  schematically  in  F'igure  1, 
where  iwt)  hot  species  of  dissimilar  mass  are  injected  into  a colder  background.  For  a 
light  species  (e  g electrons),  colliding  with  a colder,  heavy  species  (eg  hydrogen),  the  dom- 
inant process  IS  pitch-angle  scattering  Conversely  when  a heavy  species  (eg.  uranium) 
cools  due  to  collisions  with  a lighter  species,  the  dominant  process  is  loss  of  energy.  F'ig- 
ure I shows  how  these  two  cases  result  in  different  time  evolutions  for  /. 

I wo  examples  of  another  physical  phenomenon  are  depicted  in  F'igure  2.  Here 
velocity-space  is  divided  into  .1  regions,  labeled  fj,  ( and  I Region  F contains  escaping 
particles  (as  in  a mirror  machine  or  pinch).  C contains  confined  particles,  and  T contains 
trapped  particles  (as  in  a tokamak)  Boundaries  such  as  these  are  present  in  the  velocity- 
space  of  many  systems  of  interest  m reactor  studies  Such  boundaries  are  also  usually  in 
different  locations  in  F||  space  for  each  species 
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1 — V cloviu  ^^>.Kc  liinusion  o(  lighi  (1.  — cicciron)  aiul  heavy  (U  — uranium)  ions  under  Ihc  influence  ot  colli- 
smns  vMih  a void  h.ivk^tround  (H  • hydrogen)  1 he  conceninv  conlours  represeni  ^ucc'esslvc  stages  of  evolution  of 
numher  denMty  contours 


( ijs  2 — VcUKily-sp.u.c  regions  and  boundaries  The  figure  shows  typical  escape  regions  (HI  and  loroidal  trapping  (M 
regions,  with  ihc  remainder  (Kcupied  by  confined  (C)  mailer 
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Hnally.  observe  that  liquation  2 is  in  conservative  form.  Energy  losses,  such  as  radia- 
tion, appear  as  "convection''  terms  (in  1,),  since  they  conserve  the  number  of  particles.  On 
the  other  hand,  particles  may  be  added  or  removed  by  nuclear  reactions,  or  injection  or  loss 
processes  Such  features  are  treated  by  including  source  and  sink  terms  in  the  equation. 

In  brief,  this  is  the  physics  the  Fokker-Planck  Equation  can  be  used  to  study  It 
concentrates  on  the  velocity-space  dynamics  of  the  various  species  in  a plasma,  it  follows 
the  evolution  of  the  distribution  function  of  each  species.  It  can  represent  the  loss  cones 
of  mirror  machines  and  the  trapped  particles  of  tokamaks.  Injection,  losses,  and  nuclear 
reactions  can  be  included.  How  these  things  are  done  is  described  in  the  following  sections 

MATHEMATItS 


The  evolution  of  the  distribution  function  under  the  influence  of  collisions  is  a 
.Markov  process  Ihat  is,  it  is  a random  process  in  which  each  step  depends  only  on  the 
result  of  the  immediately  preceding  step.  Ihus  / is  a random  variable  of  the  coordinates 
tV,t)  which  satisfies  the  evolution  equation  ’ 

H\\i)  =//  (V  - AV.  r - At) /)(V  - A V,r;  AV  )(/■''  AV  (3) 

where  /ms  a probability  distribution  satisfy  ing 

J p (V,/;AV)  AV  = 1 <4) 


for  all  (V.t). 


fo  make  use  of  Equation  3.  we  expand  it  in  a Taylor  series, 
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where  we  have  taken  A/ and  A I sullicicntly  small  that  we  can  neglect  higher  order  terms  in 
the  expansion 


.Since  / anil  the  ditt'erential  operators  are  independent  of  V,  the  integrals  can  be  carried 
out.  yielding 
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where  l.q  (4)  has  been  used,  and  P and  Q remain  to  be  determined. 
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In  an  elegant  but  tedious  calculation,  the  properties  of  Coulomb  collisions  are  used  to 
obtain  P and  Q.'  It  is  important  to  note  that  this  calculation  assumes  both  a minimum  and  a 
ntaximum  range  for  the  Coulomb  interaction,  to  obtain  the  frequently  used  'Coulomb  loga- 
rithm" In  \ . ' 


Before  writing  down  general  expressions  lor  P and  Q,  it  is  useful  to  recall  that  the  sys- 
tems we  are  concerned  with  are  multispecics  plasmas  It  is  thus  appropriate  to  rewrite  Eq 
(6)  in  a form  which  includes  the  collisions  of  each  species  with  all  others,  including  selt- 
collisions.  It  can  be  written 


<11. , _ I 

\ii  “ '7  2 dv  d'v  dv  dv  dv  " dv 


t7) 


where 


* 

It 


5 


(8) 


.inii 


,U„(V)  = / (V)  I V - V"  I r/J  V', 

1)1  + )))f^  f 

(V  ) = - ■'  — j (V)  / 1 V - V'  I V', 

= ^ ^Ij'l  '^1 

1 2 kt2 

4t7  .V/,f 


(9) 

(10) 


In  22  ~ In  («'  ‘ 


7^ 


< 10  -“'  "A'  < /;, 


(ID 


I he  "Rosenblulh  potentials"  ,e  and  h are  integrals  over  the  distribution  functions  which 
in  principle  (see  next  section)  need  to  be  calculated  once  for  each  species,  then  summed  in 
\arious  ways  to  determine  the  effects  of  collisions. 


The  basic  Kokker-Planck  equation  (7)  is  parabolic.  Thus  Dirichlet  or  Neumann  condi- 
tions on  a closed  boundary,  plus  initial  conditions  within  the  boundary,  constitute  a well- 
poseil  n.  This  is  what  is  normally  solved  in  cases  of  physical  interest. 

•gral  operators  in  ,e  and  h cause  sufficient  difficulty  in  dealing  with  these 
c analytic  solutions  of  the  system  (7-10)  have  been  attempted  only  for  special 

^ such  case  assumes  each  species  is  Maxwellian,  with  boundary  conditions  / = 0 

at  1 = ■ .'\  more  realistic  model  is  ohiained  by  expanding  the  distribution  function  and 

integrals  m Legendre  polynomials  |l  qs.  (18-20))  and  retaining  only  the  n = 0 and  n = 1 
terms  I his  model  has  yielded  important  advances  in  the  transport  properties  of  plasmas.'’ ' 
In  less  special  cases,  numerical  assistance  is  usually  required. 

M MKRKS 


there  are  two  basically  diflerent  numerical  algorithms  needed  to  solve  the  Lokker- 
Planck  I quation  One  is  the  relatively  conventional  solution  of  the  diffusion  equation  |Lq. 
r^)|  The  other  is  the  calculation  of  the  collision  integrals  iLqs.  (8,9)1.  A cursory  examina- 
tion ol  l.qs  (7-9)  in  linite-difference  form  suggesis  that  most  of  the  work  is  involved  in 
calculating  the  integrals  l et  us  see  how  that  work  can  be  reduced.'’ 


first.  It  IS  not  necessary  to  calculate  both  integrals,  since  they  are  related  by 

JL  JL  ^ _ 

fiV  ' dV  )„~  +'n,^ 


(12) 


Thus,  once  g is  calculated.  )i  can  be  obtained  by  differencing  it.  Alternatively,  if  h is  calculat- 
ed. .g  is  the  solution  of  a Poisson  equation.  In  fact,  if  a method  for  solving  a Poisson  equation 
IS  readily  available,  it  can  also  be  used  to  obtain  //,  since  the  Laplacian  of  Kq.  (9)  can  be 
» ritten 


il\  dV 
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Belore  choosing  between  computing  the  integrals  and  solving  the  Poisson  equations,  it  is 
necessary  to  select  the  coordinate  system  to  be  used  in  velocity  space. 


Ihe  required  coordinate  system  creates  much  ol  the  trouble.  T he  spherical  sym- 
metry ol  the  C oulomb  interaction,  and  of  the  most  common  solutions  (Maxwellian),  strongly 
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suggests  the  use  of  spherical  coordinates.  Most  of  the  problems  of  physical  interest  have 
axial  symmetry,  with  the  magnetic  and  electric  fields  (if  any)  in  the  axial  direction  Note 
that  most  of  the  complexity  of  what  follows  is  due  to  this  choice  of  spherical  coordinates 

l et  us  represent  our  velocity  space  in  a spherical  coordinate  system  with 

A f 

metric  <//^  = tin  + r'dn^  + sin </</>  ^ . We  shall  assume  that  — =0.  Then  indivi- 

h<l> 

dual  elements  of  Hr.n)  represent  toroids  in  velocity  space,  and  the  integrals  (8,9)  can  be  ex- 
pressed in  terms  of  complete  elliptic  integrals 

( f ,«  ) = J"  1/ A ( iv)  /^  ( f.  H ')  ( f ")  ^ sin  W ' r/T' (14) 

and 

III  + III.  p 

=~ I //,()".  W)  (r)^  sin  «'r/r  r/H'.  (15) 

tin  j n 

where 

u = 4 | 2 + (^  ■)2  _2  f'T  (cos  « cos  ft'  - sin  ft  sin  «') ‘'^  (16) 

W =16  f f "sin  « sin  ft'  / fy'  (17) 

If  / is  represented  on  an  /V,  x N„  grid,  then  the  kernel  of  these  integrals,  e.g. 

ft;  r,  ft  ) = L'A'(  ff)  r sin  ft', 

is  an  X table  of  values  for  each  pair  of  species,  requiring  an  ,V,  x A',  dot  product  for 
each  point  and  pair  of  species. 

There  is  an  alternative  which  appears  to  reduce  the  work  involved  in  computing  these 
integrals.^’’  It  is  based  on  methods  which  have  been  developed  for  the  Poisson  equation  Sup- 
pose we  expand  the  distribution  function  in  Legendre  polynomials. 


/,(L.  ft)  = X >h  (cos  ft)  (18) 
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Ihen  the  radial  and  a/imuthal 

integrals  are 

separable,  and 

the  a/imuthal  integrals  have  a 

direct  quadrature,  giving 

4t7 

f / III  1(2  . 

f . 

1 

1 

"-y  L.  2 

(19) 

4ir  -1 

J 1 1 t ) 

n 

T ~ 

1 

\ (y  .ii) 

(L)^ 

(/r 

i>, 1,(1. II) 

An 

f (...■> 

n 

r. 

/f, 

/ 1-"  \ ^ t " \ * r/  lx  ' 

(20) 

* 2/;  + 1 

1 r ,/M  1 r ; (/y 

where 

1 = min  ( L, 

f ) 

(21) 

f . = max  ( T,  f')  (22) 

In  this  way,  the  calculation  of  n and  h has  been  reduced  from  complete  x integrals 
to  one  ,V,  integral  for  each  a/imuthal  mode,  and  the  kernels  of  the  integrals  are  (frac- 
tional) polynomials,  rather  than  elliptic  functions.  Incidentally,  the  ii=0  equations  for  k 
and  li  are  the  one-dimensional  (radial)  integrals  for  the  Fokker-Planck  Lquation. 
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I he  I.egcmJre  expansion  is  a greai  improvement  if  only  a few  a/imuthal  modes  are 
needed  to  describe  the  distribution  of  interest.  However,  il  localized  distributions  occur,  for 
example  when  intense  beam  injectors  are  present,  a number  of  modes  may  be  re- 

guired  Then  the  work  can  be  comparable  to  direct  integration  (elliptic  integrals  are  easy 
to  calculate) 

There  is  no  consensus  on  this  choice.  Some  workers  solve  Tigs.  (14-17)  while  others 
solve  Tgs  (14-22)  The  choice  is  to  some  extent  dependent  on  the  physical  problem  being 
.iddressed  At  the  present  lime,  the  direct  solution  appears  to  have  an  advantage  for  computa- 
tional reasons  It  is  a very  highly  ordered  operation  on  a lot  of  data,  and  as  such  its  computa- 
tional cost  IS  reduced  up  to  two  orders  of  magnitude  on  a vector  computer.  (.See  lable  I.)  In 
any  case,  calculating  these  integrals  is  only  half  o(  the  problem. 

Ihe  other  part  ol  the  Tokker-Planck  quadrature  is  solving  the  diffusion  equation  (Tiq. 
(7l|  In  spherical  coordinates,  expansion  of  the  operations  in  Tiq.  (7)  is  tedious,  but 
straightforward  II  I q (7)  is  transidrmed  into  the  Legendre  representation  of  /,  the 
ditferential  operators  in  n are  replaced  by  coupling  terms  between  the  various  modes. 

Since  both  /and  the  integrals  g and  /i  have  azimuthal  dependences,  the  spectrum  ol 
the  operator  on  the  right-hand  side  of  Tq.  (7)  is  not  closed;  that  is,  if  / contains  modes  only 
up  to  ,V,  Tquation  (7)  can  yield  modes  in  with  n > fj.  Some  computational 

method  must  then  be  prescribed  to  truncate  the  spectrum  of  the  operators.  That  is,  the 
part  of  / which  is  "lost”  to  higher  modes  must  be  returned  to  Ihe  lower  modes,  in  a way 
which  conserves  mass,  etc.  as  far  as  possible.  This  is  another  problem  with  the  l.egendre 
representation  Because  of  these  dilTiculties.  the  following  remarks  will  be  restricted  to 
the  dillerence  formulation  of  Lq.  (7). 

T.quation  (7)  contains  both  a convective  and  a dilTusive  part,  so  several  numerical  sta- 
bility criteria  must  be  satisfied  by  any  explicit  algorithm  for  solving  it.  In  terms  of  the  (7  and  r 
ol  Tq  (6),  they  are  approximately 

< y „„„  ,2„ 


l>  1 


p US  some  obscure  and  nasty  ones  on  the  oH-diagonal  terms.  The  usual  case  is  that  the 
l.isi  term.  • ‘nti ‘’//i,,,,  is  the  most  restrictive,  in  the  cells  near  the  origin  of  the  coordinate 
system  As  is  typical  in  computational  lluid  dynamics,  the  stability  considerations  for  an 
algorithm  .ire  dominated  by  a region  in  which  very  little  interesting  physics  is  going  on.  In 
most  vases  of  physical  interest,  the  <M  required  for  explicit  stability  is  too  small  to  be  prac- 
iic.il,  ,md  an  implicit  algorithm  must  be  used 

Il  IS  clear  Irom  these  considerations  that  the  grid  should  not  contain  a cell  at  f =0. 
Ihe  most  natural  choice  of  linite-dillerence  cells  uses  the  origin  and  the  axis  as  cell 
boundaries  I hen  the  axial  boundary  conditions  take  care  of  themselves:  no  llux  Hows 
through  these  boundaries 

W hen  the  ilillerential  operators  in  Lq.  (7)  are  written  out  in  components,  a lot  of  terms 
appear  ‘ Keduvtion  to  dillerence  form  requires  a lot  of  algebra,  but  is  otherwise  straightfor- 

w.ird  Ihe  speed  ol  present-day  computers  makes  practical  what  was  prohibitive  a few 


I 


years  ago  Kcpiaeing  ihe  operator  on  the  right-hand  side  of  liq.  (7)  by  centered  differences 
results  m a ‘J-point  spatial  difference  operator  (5-pomt  if  second-order  accuracy  is  not  re- 
tained) 

X I ni  + I.  I + m)  = Riij)  (25) 

I - I II. I 
I II  I 

■\d\ancing  / in  time  then  requires  solving  a (,V,  x ,V„  )’  sparse  matrix  equation  Typical  griil 
resolution  is  .V,  = 50,  .\„  =20,  so  in  principle  a 1000x1000  matrix  inversion  is  involved 
In  lad,  methods  are  novx  available  for  solving  such  banded  n.airix  equations  by  rapidly  con- 
\ergent  iterations,**,  or  cyclic  reduction.'^ 

lahle  I illustrates  the  gams  available  with  such  methods.  It  presents  speeds  obtained  on 
the  \Rl  (2-pipe)  Texas  Instruments  ".Advanced  Scientific  C omputer."  The  two  most  expensive 
elements  of  the  I okker-l’lanck  solution  are  the  integral  of  f-.q.  (8)  and  Ihe  solution  of  the 
matrix  equation  It  q (2,‘')|  Table  I presents  the  scalar  time,  vector  time,  the  ratio  of  these 
limes,  and  the  vector  speed  in  ntillions  of  lloating-pomt  operations  per  -cond  (MFT.OPS)  for 
these  tasks  To  put  this  speed  in  perspective,  performing  the  integrals  of  l:q  (8)  for  each 
point  ol  a 20x511  grid  requires  2 million  operations,  anil  lakes  40  n. -.'.liseconds,  about  the 
lime  It  lakes  to  write  one  word  on  a disk  (disk  rotation  lime) 

Initial  and  lexterior)  boundary  conditions  are  usually  fairly  simple.  T he  initial  condi- 
tions usually  rapidly  disappear  behind  the  dynamics  ol  the  diffusion  There  are  no  wave- 
like  solutions  for  them  to  excite  They  are  typically  Vfaxwellian  distributions,  possibly 
with  some  delta-liinction  beams. 

Boundary  conditions  require  more  care.  .At  large  velocities,  collisions  become 
less  frequent,  so  particles  which  get  there,  can  stay  there  In  particular,  if  an  electric  tield  is 
present,  there  is  a velocity  (lor  each  species)  above  which  a particle  "runs  away."  I his 
means  that  the  final  slate  ol  such  a calculation  is  one  in  which  all  particles  have  run  away  , 
but  lor  practical  problems,  the  lime  it  takes  this  to  happen  is  long  compared  to  experimen- 
tal limes  Specific  values  may  be  prescribed  at  the  (exterior)  boundary,  or  it  may  be  taken 
impermeable  Whatever  choice  is  made,  the  parameters  of  the  calculation  can  usually  be 
chosen  so  the  results  of  interest  arc  insensitive  to  them 

finally,  it  should  be  mentioned  that  there  are  many  "collisional"  processes  that  are 
not  simple  ( oulomb  processes.  T or  example,  the  particles  may  be  scattered  by  electros- 
tatic plasma  waves  Such  coherent  processes  are  not  described  by  Ihe  Coulomb  integrals 
II  qs  (8.9)1,  but  the  diffusion  and  drag  resulting  from  them  can  readily  be  included  in  the 
I okker-IManck  formalism 

AIMM.KMTONS 

Slight  modifications  of  the  f okker-l’lanck  liquation  are  often  needed  to  apply  it  to 
specific  experimental  configurations.  This  section  indicates  what  changes  are  required  to  study 
typical  mirror,  pinch,  tokamak.  or  laser-plasma  problems.  These  applications  are  il- 
lustrative of  the  adaptability  of  velocity-space  methods,  and  show  how  details  of  an  ex- 
periment can  be  represented. 

The  dominant  feature  of  a mirror  machine  is  the  presence  of  loss  cones  in  velocity 
space  *■  I his  means  that  particles  m a region  such  as  li  of  ITgure  1 leave  the  machine  in  a 
very  short  time  Thus  the  boundary  condition  1 — 0 should  be  applied  on  the  surface 
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I.iblc  1 

Table  I presenis  measured  perlbrnianee  data  lor  ibe  most  expensive  computational  elements  of 
the  1 okker-l’lanek  calculation  I quaiion  (8)  is  readily  veclori/able,  while  Hq  (25),  which  in- 
\ol\es  recursive  calculations,  is  not  I he  "scalar"  and  ''vector”  speeds  were  measured  on  a 
2-pipe  Texas  Instruments  AS(  , and  are  expressed  in  terms  of  "millions  of  floating-point  opera- 
tions per  secoiul  ' 
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between  regions  C'  and  f . What  is  unlike  other  lluid  problems  is  that  the  boundary  sur- 
lace  ( -I  IS  dillerent  for  each  q/m.  for  example,  a "square  well"  magnetic  mirror  with  mir- 
ror ratio  and  maximum  potential  'l>  has  loss  cones 


sin  ' It  = 


2^1' 

nil^ 


(2b) 


which  are  dillerent  for  each  particle  q,  m. 


fhe  methods  used  to  solve  these  boundary  value  problems  are  those  of  conventional 
tluid  mechanics,  there  arc  just  more  boundaries  than  usual.  The  boundaries  may  be  treated 
properly,  solv  mg  f q (7)  [that  is,  (!q.  (25)1  subject  to  them  quick  and  dirty  alternative  is  to 
solve  the  problem  with  boundary  conditions  at  ) = oo  (or  some  ),  and  then  set  / =0  in 
region  f at  each  time  step,  for  a dilfusive  problem,  this  will  not  result  in  an  instability, 
just  a lirst-order  error  in  nr;  but  if  an  implicit  method  is  in  use  which  is  first  order  in  nr, 
such  terms  are  already  present.  Typical  computational  results  include  the  loss-rate  of 
ions,  and  the  potential  ■!>  necessary  to  obtain  charge  neutrality. 


pinch  experiment  is  characterized  by  a rapid  change  of  the  magnetic  Held  B in 
time  fhe  appropriate  t okker-Planck  liquation  is  obtained  by  replacing  the  left-hand  side 
ol  I q (71  by  the  left-hand  side  o(  I q (2)  The  change  with  time  of  /(produces  a convection 
of  / in  the  perpendicular  direction  (Kecall  that  W is  aligned  along  the  axis  of  our  coordinate 
system  ) A lluid  element  convects  at  a rate 


</l  ^ 

ill 


1 A 

2 H ill 


(27) 


The  (actor  of  1/2  is  present  because  there  are  two  components  of  Tin  the  perpendicular 
direction  The  physical  (irocess  involved  is  constancy  of  the  magnetic  moment.  In  typical 
pinch  experiments,  this  compression  heats  the  ions.  The  f'okkcr-Planck  liquation  is  used  to 
obtain  the  resulting  shape  of  /.  This  determines  the  hot  ion  lifetimes  and  their  rate  of  isotrop- 
i/ation,  i.e.,  the  rate  at  which  the  ions  acquire  parallel  velocity  and  are  lost  from  the  pinch. 
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I okamaks  ha\e  lour  I'eatures  v^-hith  arc  important  to  a vclocity-spacc  model.  First, 
an  electric  field  is  present,  parallel  to  the  magnetic  field,  the  plasma  carries  a net  current 
Second,  there  are  magnetic  mirrors  within  each  magnetic  surface,  and  they  can  trap  some 
particles.  Third,  there  is  i significant  variation  of  plasma  conditions  with  spatial  posi- 
tion: velocitv -space  modelling  of  a tokamak  is  viable  because  individual  magnetic  surfaces  are 
fairlv  well  isolated  from  each  other,  and  conditions  on  each  magnetic  surface  are  fairly  uni- 
form f ourth.  there  is  a relatively  high  level  of  impurities  in  most  tokamaks.  N'arious  work- 
ers have  examined  one  or  more  of  these  features.  Thev  have  not  yet  all  been  treated  together 

I he  electric  field  results  in  distortion  of  the  distribution  function,  and  rapid 
convection  along  the  axis.  Ihis  can  produce  difficulties  for  the  algorithm  which  solves  fq 
(25).  particularly  around  the  origin.  Physical  problems  of  interest  are  the  distortion  of  / and 


the  resultant  resistivity:  the  ratio  of  the  current  it  carries  to  E. 


The  calculation  of  the 


resistivitv  has  been  explored  thoroughly  when  no  toroidal  effects  are  present  Present  work 
thus  involves  inclusion  of  some  of  the  other  effects  as  well 

Trapped  particle  effects  can  be  included  by  treating  separately  the  trapped  and  un- 
trupped  particles.  (Regions  1 and  C respectively  of  figure  I)  The  particles  in  C satis- 
fy a conventional  f okker-Planck  Eguation,  while  those  in  T satisf'  kokker-Planck  equa- 
tion with  no  electric  field,  since  they  are  not  free  to  accelerate  in  that  iield.  and  the  cr.’.Tgy 
gam  in  one  bounce  is  negligible  compared  to  their  average  energy.  This  can  be  most  sin'niy 
modeled  by  making  /:  a function  of  T.-vehich  is  the  lokamak  value  in  C.  and  /era  in  T Cal- 
culations with  this  model  yield  the  toroidal  corrections  to  the  resistivity.'^ 

I he  magnetic  surfaces  are  really  not  independent,  since  ion  orbits  (eg.,  banana  or- 
bits) extend  over  different  surfaces.  Convection  and  loss  mechanisms  also  couple  surfaces 
Solutions  of  a coupled  radial  transport  equation  plus  a f okker-Planck  equation  on  several  sur- 
laces  have  been  atlempled.'*  however  a large  number  of  approximations  are  required.  This 
IS  siill  a very  dilficult  problem 

f inally,  impurities  introduce  two  important  effects,  first,  they  greatly  ’increase 
pitch-angle  scattering  ol  trapped  particles,  beams,  etc.  .Second,  if  they  have  high  atomic 
numbers,  they  radiate,  resulting  in  an  electron  energy  loss.  The  first  effect  is  important 
when  studying  beam  injection  The  second  is  important  in  overall  energy  balance 

As  a last  applicaiion,  Ihe  f okker-Planck  fqualion  has  several  interesting  areas  ol 
ihe  laser-plasma  interaction  to  investigate  In  Ihe  absorption  region,  where  the  laser  elec- 
tric field  IS  largest,  Ihe  electron  distribution  function  is  swept  back  and  forth  in  veUKily  space 
I rom  the  point  ol  view  of  the  ions,  the  electron  distribution  function  acquires  an  ellipti- 
cal dislorlion  f he  fokker-Planck  f qualion  can  be  used  to  quantitatively  calculate  the 
healing  rale  due  to  ihis  process 


Both  in  Ihe  lokamak  and  laser  plasma  environments,  beams  are  formed  in  one  region 
and  transported  to  (iniected  into)  another  The  f-okker-Planck  [.qualion  can  calculate  Ihe  in- 
ler.iciion  between  these  beams  and  Ihe  background  through  which  they  are  propagating 
The  beams  are  simply  sources  in  the  equation. 


i| 


, » 
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I hose  .ipplioations  indioaie  both  the  \ersatiliiv  ol  the  f okker-Planek  equation,  and  the 
amount  i.!  elTort  vshieh  is  required  to  solve  such  problems  The  computational  cost  of  these 
solutions  IS  rapidK  decreasint:.  and  more  examples  of  such  applications  are  appearing  Will 
>ours  be  the  next  ’ 

RK  \(  TORS 

I’  ma\  .ippear  Irom  the  previous  section  that  nuclear  reactions  have  been  neglected 
m this  \elocii\ -space  work  That  is  not  so  The  applications  section  showed  how  velocity- 
space  methods  are  adaptable  to  the  physical  characteristics  ol'  various  devices  In  most  of 
the  experiments  describeil  above,  reactor  studies  have  been  carried  out  with  the  resulting 
model  Ihe  nuclear  reaction  calculation  is  usually  the  most  conventional  part  ol  such  calcu- 

l.iiions  It  uses  the  calculated  distribution  runctions  to  compute  reaction  rales  from 
the  .ippropri.ite  cross-seclioiis 

Once  the  healing  ,ind  loss  processes  of  the  experiment  are  modeled,  it  is  easy  to  inject 
the  reading  species,  aiul  c.ilculaie  their  dynamics.  One  part  of  that  dynamics  is  the  nu- 
clear reactions  they  undergo  The  reaction  integral  determines  two  things;  the  rale  of  pro- 
diKtion  ol  nucle.ir  energy,  and  the  reaction  conslitueni  sources  and  sinks  in  velocity  space 

\t  this  point,  you  m.iy  take  you,’’  choice  and  pay  your  money  There  is  no  argu- 
ment about  the  neutrons  I hey  give  their  energy  to  the  external  environment  The  fate  of 
the  reaction  products  is  sonreihing  else  Ihe  inexpensive  choice  is  to  deposit  the  ene'^gy  of 
the  charged  reaction  priidiicis  where  it  will  probably  go  first;  namely  , into  the  electrons 

Ihe  expensive  choice  is  to  lollow  the  velocity-space  behavior  of  these  products.  It 
IS  expensive  lor  two  reasons  First,  the  velocities  with  which  the  reaction  products  are  born 
IS  very  l.irge  compared  to  the  original  plasma  velocities.  Thus  the  required  velocity- 
coonlinaic  -.pace  becomes  larger.  ,md  liner  griddmg  is  required.  Second,  collision  in- 
tegrals are  required  lor  each  p.iir  ol  species,  so  much  more  work  is  needed  for  each  addition- 
al reaction  product  which  is  treated 

H'ji  It  IS  the  expensive  choice  that  will  be  needed  in  the  coming  years  In  practice,  all 
the  proposed  plasma  re.ictors  h.ive  loss  regions,  and  it  is  important  to  know  how  much  en- 
ergy IS  returned  to  the  reactor  (dasma  before  products  arc  cither  lost,  or  slow  down  This  will 
also  .illow  a quaMiilative  .issessment  o)  nuclear  reactions  in  which  the  participants 
themselves  participate  Such  calculations  are  nowon  ihe  drawing-boards,  wailing  their  turn 

This  IS  presently  a wide-open  area  ol  research,  and  even  a factor  of  two  improvement 
in  yield  lor  some  experimental  design  would  be  important  news  Designs  such  as  the 
■ iw o-component  torus"  are  presently  on  the  edge  of  scientific  breakeven,  and  looking  for 
any  available  help 

This  survey  has  indicated  what  kinds  ol  problems  can  be  addressed  by  velocity-space 
methods,  and  how  they  are  modeled  The  basic  conclusion  is  that  all  of  the  major  plas- 
ma reactor  designs  are  amen.ibic  to  such  treatment,  and  each  has  had  some  important  as- 
pect analyzed  in  this  way  It  is  .ilso  clear  ihai  much  remains  to  be  done  The  cost  of 
velocity -space  computaiions  has  been  large  m the  past,  but  is  now  coming  down.  Scientists 
j ; who  have  an  acquaintance  with  distrihution-lunction  methods  should  find  a friend  m the 

i I okker-l’lanck  I quation 
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